/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains the class blackbody_grain_model.

// $Id$

#ifndef __blackbody_grain__
#define __blackbody_grain__

#include "mcrx-types.h"

#include "constants.h"
#include "mono_poly_abstract.h"
#include "grain_model.h"
#include "blackbody.h"

namespace mcrx {
  template <template<class> class, typename> class blackbody_grain_model;
  template <template<class> class chromatic_policy, typename T_rng_policy> 
  class blackbody_grain_model_creator_functor;
};

/** This grain model consists of blackbodies. There is only one
    parameter and that is the opacity, which is independent of
    wavelength. If you assume that the grains are spheres of a certain
    radius and density, you discover that the opacity is given by
    \f[\kappa = \frac{3}{4\rho R},\f] which enables you to compute the
    emitted luminosity per mass. You get \f[L=4\kappa\sigma_{SB}T^4
    M.\f]

    To get the output SED, we integrate the blackbody intensity over
    area and solid angle, and we get \f[L_\lambda/M = A \pi B_\lambda
    /M = \frac{4 \pi^2 R^2 B_\lambda}{4\pi R^3/3} = 4\pi \kappa
    B_\lambda.\f] Integrating this over wavelength returns the
    expression above.

    The determination of the equilibrium temperature is similar to
    what's done for thermal_equilibrium_grain, except here it's much
    simpler. The equilibrium temperature does not depend on the
    assumed radius since the area cancels when calculating the
    absorbed and emitted luminosities when the cross-section is not
    wavelength-dependent. For one grain, \f[L_\mathrm{abs} = 4 \pi
    \int I_\lambda \sigma d\lambda = 4 \pi \pi R^2 \int I_\lambda
    d\lambda = L_\mathrm{emitted} = 4 \pi R^2 \sigma_{SB} T^4,\f] and
    you end up with \f[T^4 =\frac{\pi \int I_\lambda
    d\lambda}{\sigma_{SB}}.\f] See the documentation for
    thermal_equilibrium_grain for more info.
*/
template <template<class> class chromatic_policy, 
	  typename T_rng_policy = mcrx::global_random> 
class mcrx::blackbody_grain_model : 
  public mcrx::grain_model<chromatic_policy, T_rng_policy> {

private:
  const T_float kappa_; ///< The opacity (cross section/mass) of the grains.
  array_1 alambda_; ///< The wavelength vector for absorption calculations.
  array_1 elambda_; ///< The wavelength vector for emission calculations.

public:
  blackbody_grain_model(T_float k, T_unit_map u) : kappa_(k) {
    // copy units from those supplied into *scatterer* unit map.
    T_unit_map& u1 = 
      mcrx::scatterer<chromatic_policy, T_rng_policy>::units_;
    u1["length"] = u.get("length");
    u1["wavelength"] = u.get("wavelength");
    u1["mass"] = u.get("mass");
    T_unit_map& u2 = 
      mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
    u2["length"]=u.get("length");
    u2["wavelength"]=u.get("wavelength");
    u2["mass"] = u.get("mass");
  };
  blackbody_grain_model(const blackbody_grain_model& rhs) :
    kappa_(rhs.kappa_), alambda_(independent_copy(rhs.alambda_)),
    elambda_(independent_copy(rhs.elambda_)) {}

  virtual boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> > 
  clone() const {
    return boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> >
      (new blackbody_grain_model(*this));};

  const array_1& alambda() const {return alambda_;};
  const array_1& elambda() const {return elambda_;};

  /** Resamples all the grain objects to use the specified wavelength
      vector. Since there is no wavelength info, it's pretty simple. */
  void resample (const array_1& lambda, const array_1& elambda=array_1()) {
    const array_1 el(elambda.size()>0?elambda:lambda);
    assign(alambda_, lambda); 
    assign(elambda_, el);
    array_1 opacity(lambda.size());
    opacity=kappa_;
    // assign to the scatterer members. albedo and g are set to 0.
    assign(this->opacity_, opacity);
    assign(this->albedo_, 0.0*opacity);
    this->phfunc.set_g (array_1(0.0*opacity));
  };

  /** This redefinition of the virtual scatterer::set_wavelength
      function just calls the resample function (which also
      recalculates the extinction curve). */
  virtual void set_wavelength (const array_1& lambda) {
    resample(lambda);
  };

  /** This gets called by the template virtual hack in
      grain_model::calculate_SED. */
  template <typename T>
  void calculate_SED_virtual(const blitz::ETBase<T>& intensity, 
			     const array_1& m_dust, array_2& sed,
			     array_2* temp) const;

  /** Calculates the temperature corresponding to the specified
      intensity vector. */
  template <typename T>
  T_float calculate_temp(const blitz::ETBase<T>& intensity) const;
};


template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
template <typename T>
void
mcrx::blackbody_grain_model<chromatic_policy, T_rng_policy>::
calculate_SED_virtual(const blitz::ETBase<T>& intensity, 
		      const array_1& m_dust,
		      array_2& sed, array_2* temp) const
{
  using namespace blitz;

  if(temp)
    // temp does not make sense for this grain
    *temp = blitz::quiet_NaN(T_float());

  const T& i(intensity.unwrap());
  const int nc=m_dust.size();
  const array_1 invelambda(1./elambda_);
  array_1 tempsed(elambda_.shape());

  // convert kappa*m_dust to SI units (mass unit cancels, only area left)
  const T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  const T_float k = 
    kappa_*units::convert(u2.get("length")+"*"+u2.get("length"),"m^2");
  // gdb can't see variables defined within block... pos
  T_float i_int, temp4, B_int, scalefactor;
  for(int c=0; c<nc; ++c) {
    if(m_dust(c)==0)
      continue;

    // this is so much simpler than for a wavelength-dependent opacity...
    i_int =
      integrate_quantity(i(c,Range::all()), alambda_, false);

    temp4 =
      constants::pi*i_int/constants::sigma_SB;

    // make sure we preserve the luminosity as realized on the grid
    // (this is not a good idea if the wavelengths don't cover the
    // full emission, but is by necessity how the thermal equilibrium
    // grain works)
    tempsed = 4*constants::pi*k*m_dust(c)*B_lambda(invelambda, pow(temp4,-0.25));

    B_int = 
      integrate_quantity(tempsed, elambda_, false);

    scalefactor = 4*k*constants::sigma_SB*temp4*m_dust(c)/B_int;
    DEBUG(2,cout << "Rescaling grain luminosity by " << scalefactor << endl;);

    // can't scale NaNs or infs, they both mean vanishingly small lum
    if ((scalefactor==scalefactor)&&(scalefactor<blitz::huge(T_float()))) {
      sed(c,Range::all())+= scalefactor*tempsed;
    }
    ASSERT_ALL(sed(c,Range::all())>=0);
  }
}


template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
template <typename T>
mcrx::T_float
mcrx::blackbody_grain_model<chromatic_policy, T_rng_policy>::
calculate_temp(const blitz::ETBase<T>& intensity) const
{
  const T& i(intensity.unwrap());

  // convert kappa*m_dust to SI units (mass unit cancels, only area left)
  const T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  const T_float k = 
    kappa_*units::convert(u2.get("length")+"*"+u2.get("length"),"m^2");

  const T_float i_int =
      integrate_quantity(i, alambda_, false);
  const T_float temp =
    pow(constants::pi*i_int/constants::sigma_SB, 0.25);
  return temp;
}


template <template<class> class chromatic_policy, typename T_rng_policy> 
class mcrx::blackbody_grain_model_creator_functor
{
public:
  typedef boost::shared_ptr<mcrx::grain_model<chromatic_policy, T_rng_policy> > return_type;
  
  blackbody_grain_model_creator_functor() {};
  return_type operator()(const Preferences& p) const {

    //THIS IS A HACK BECAUSE WE CAN ONLY HAVE ONE ARGUMENT
    T_unit_map units;
    units["mass"] = p.getValue("massunit", std::string());
    units["length"] = p.getValue("lengthunit", std::string());
    units["wavelength"] = p.getValue("wavelengthunit", std::string());

    const T_float k = p.getValue("blackbody_opacity", T_float());
    
    return return_type
      ( new mcrx::blackbody_grain_model<chromatic_policy, T_rng_policy> 
	(k, units));
  };
};


#endif

